########################################################
## Replication of Carrubba et al (2008)
## Iyer, Malekovic, Solano
########################################################
##First part of code: replication; second part: our own analyses

########################################################
## Packages
########################################################

# Load 'foreign' package, so we can read STATA files
library(foreign)

# Load Zelig
library(Zelig)

# Other packages needed
library(sandwich)
library(lmtest)


########################################################
## Loading the Data
########################################################

data <- na.omit(read.dta("apsrreplication.dta"))

# Note: for some reason, not all NAs are invalid
data.isNotNA <- is.na(data[,"ECJPlAgree"]) & (data[,"normnetwobs"] != 0)
data[data.isNotNA,"ECJPlAgree"] <- 0
data <- na.omit(dataNA)

# Defines default cluster, and helper method to call vcov
vars.G <- c("casenumber")

# Numeric control variables
vars.Z <- c("CommIsDef", "CommIsPl", "CommObsPl", "CommObsDef", "percham")
data.Z <- data[,vars.Z]

# Representative means
data.Z.means <- apply(X = data.Z, MARGIN=2, FUN=mean)
data.Z.medians <- apply(X = data.Z, MARGIN=2, FUN=median)

########################################################
## Custom Functions
########################################################

## Implements mode in R
my.mode <- function(x) as.numeric(names(sort(-table(x)))[1])

## Extract a robust, clustered variance-covariance matrix         
my.vcov <- function(fm, cluster){
         M <- length(unique(cluster))   
         N <- length(cluster)           
         K <- fm$rank                        
         dfc <- (M/(M-1))*((N-1)/(N-K))
         # dfc <- M / (M - 1) // Another possible finite sample adjustment    
         uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
         vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
}

## Implements clustered vcov by default
my.vcov.c <- function(fm) my.vcov(fm, fm$data[,vars.G])

## Implements robust standard errors (people.su.se/~ma/clmclx.R)
my.se <- function(fm, cluster) coeftest(fm, my.vcov(fm, cluster))

## Implements robust standard errors (people.su.se/~ma/clmclx.R)
my.se.c <- function(fm) coeftest(fm, my.vcov.c(fm))


## Extract an alternate robust, clustered variance-covariance matrix
my.vcov2 <- function(fm, cluster){
         M <- length(unique(cluster))   
         N <- length(cluster)           
         K <- fm$rank                        
         ## dfc <- (M/(M-1))*((N-1)/(N-K))
         dfc <- M / (M - 1) 
         uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
         vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
}

## Implements alternate clustered vcov by default
my.vcov2.c <- function(fm) my.vcov2(fm, fm$data[,vars.G])

## Implements alternate robust standard errors (people.su.se/~ma/clmclx.R)
my.se2 <- function(fm, cluster) coeftest(fm, my.vcov2(fm, cluster))

## Implements alternate robust standard errors (people.su.se/~ma/clmclx.R)
my.se2.c <- function(fm) coeftest(fm, my.vcov2.c(fm))


## Sets values
my.setx <- function(X, X.set, vars = colnames(X), fun=mean, funs = list(), values = list()) {
	if (is.null(X.set)) X.set <- X[1,vars]
	X.set[vars] <- apply(X = X[vars], MARGIN=2, FUN=fun)
	if (length(funs) > 0) {
		for (i in 1:length(funs)) {
			X.set[names(funs[i])] <- apply(X = X[names(funs[i])], MARGIN=2, FUN=funs[[i]])
			}	
	}
	if (length(values) > 0) X.set[,vars] <- replace(X.set[,vars], names(values), values[names(values)])
	unlist(X.set)
}

## Creates appropriate matrices of values for my.Xs from a list.
## e.g. my.list = list("normnetwobs" = c(0,0.1), "normwobsgov" = c(0,0.1))
my.values <- function(my.list) {
	do.call(cbind, my.list)
}
m22.govlit.values <- matrix(data = c(0,0.1,0,0.1), nrow=2,ncol=2)
colnames(m22.govlit.values) <- c("normnetwobs", "normwobsgov")

## Creates subsets of interest for a given representative X, and a matrix of values:
## 	 - colnames representing each variable to replace
##   - cols store the values of each variable for a given subset
##   - rows store the values across subsets
my.Xs <- function(X.set, values) {
	Xs.data <- c()
	for (j in 1:dims(values)[1]) {
		Xs.data <- c(Xs.data,as.vector(replace(X.set, colnames(values), values[j,])))
	}
	matrix(data = Xs.data, ncol=length(X.set), nrow=dims(values)[1], byrow=TRUE)
}

## Runs simulation
my.sim <- function(fm, k, Xs) {
	beta.draws <- mvrnorm(k, mu = fm$coefficients, Sigma = my.vcov.c(fm))
	pii.sim <- matrix(data = NA, ncol=dims(Xs)[1], nrow=k)
	for (j in 1:dims(Xs)[1]) {
		x <- as.vector(Xs[j,])
		pii <- c()
		for (i in 1:k) {
			pii[i] <- fm$family$linkinv(x %*% beta.draws[i,])
		}
		pii.sim[,j] <- pii
	}
	pii.sim
}

modefunction <- function(x){
	as.numeric(names(sort(-table(x)))[1])}
	
########################################################
## Descriptive Statistics - Table 1
########################################################

vars.summary <- cbind(data$ECJPlAgree, dataNA$normnetwobs, dataNA$govislit, data$article169, dataNA$article177, data$percham, dataNA$CommIsPl, dataNA$CommIsDef, dataNA$CommObsPl, dataNA$CommObsDef)

descriptive <- cbind(apply(vars.summary, 2, mean),apply(vars.summary, 2, sd), apply(vars.summary, 2, min), apply(vars.summary, 2, max))
rownames(descriptive) <- c("ECJ ruling for plaintiff", "Net weighted observations for plaintiff", "Government litigant", "Infringement Proceeding", "Article 234 case", "Chamber size", "Commission is plaintiff", "Commission is defendant", "Commission observation for plaintiff", "Commission observation for defendant")
colnames(descriptive) <- c("Mean", "Standard Deviation", "Minimum", "Maximum")

descriptive
########################################################
## Note on dummy variables
########################################################

## Rulings can be either:
## 	A- Preliminary Rulings (i.e. brought under Art. 234)
##  B- Others
##      i- Infringement proceedings (Article 226)

## The dummy variables are coded as follows:
## article177 captures whether a preliminary ruling (1) or not (0)
## article169 captures whether an infringement proc. (1) or not (0)

########################################################
## Table 2 - Model 1
########################################################

# Variables
vars.m21.Y <- c("ECJPlAgree")
vars.m21.X <- c("normnetwobs" , "percham" , "CommIsPl" , "CommIsDef" , "CommObsPl", "CommObsDef")
vars.m21 <- c(vars.m21.Y, vars.m21.X)

# Subsets the data
data.m21 <- data[,c(vars.m21, vars.G)]

# Creates formula (easier to deal with)
formula.m21 <- as.formula(paste(c(vars.m21.Y,paste(unlist(vars.m21.X), collapse="+")), collapse="~"))

# probit estimation, using glm
m21.glm <- glm(formula.m21, family=binomial(link = "probit"), data = data.m21, na.action = na.pass)
summary(m21.glm)

# calculates robust standard errors, clustered by case
m21.my.se <- my.se(m21.glm, data.m21[,vars.G])
m21.my.se


########################################################
## Table 2 - Model 2
########################################################

vars.m22.Y <- c("ECJPlAgree")
vars.m22.X <- c("normnetwobs" , "govislit", "normwobsgov", "percham" , "CommIsPl" , "CommIsDef" , "CommObsPl", "CommObsDef")
vars.m22 <- c(vars.m22.Y, vars.m22.X)

data.m22 <- data[,c(vars.m22, vars.G)]

formula.m22 <- as.formula(paste(c(vars.m22.Y,paste(unlist(vars.m22.X), collapse="+")), collapse="~"))

m22.glm <- glm(formula.m22, family=binomial(link = "probit"), data = data.m22, na.action = na.pass)
summary(m22.glm)
m22.my.se <- my.se(m22.glm, data.m22[,vars.G])
m22.my.se

#putting together models 1 and 2 into the same table
mod1 <- c(m21.my.se[2,1],m21.my.se[2,2],NA,NA,NA,NA,m21.my.se[3,1],m21.my.se[3,2],m21.my.se[4,1],m21.my.se[4,2],m21.my.se[5,1],m21.my.se[5,2],m21.my.se[6,1],m21.my.se[6,2], m21.my.se[7,1],m21.my.se[7,2],m21.my.se[1,1],m21.my.se[1,2],nrow(data.m21))

mod2 <- c(m22.my.se[2,1],m22.my.se[2,2],m22.my.se[4,1],m22.my.se[4,2],m22.my.se[3,1],m22.my.se[3,2],m22.my.se[5,1],m22.my.se[5,2],m22.my.se[6,1],m22.my.se[6,2],m22.my.se[7,1],m22.my.se[7,2],m22.my.se[8,1],m22.my.se[8,2],m22.my.se[9,1],m22.my.se[9,2],m22.my.se[1,1],m22.my.se[1,2], nrow(data.m22))

table2 <- cbind(mod1,mod2)
colnames(table2) <- c("Model 1", "Model 2")
rownames(table2) <- c("Net weighted observations for plaintiff", "se1", "Net weighted observations for plaintiff x Government litigant", "se2", "Government litigant", "se3", "Chamber size", "se4", "Commission is plaintiff", "se5", "Commission is defendant", "se6", "Commission observation for plaintiff", "se7", "Commission observation for defendant", "se8", "Constant", "se9","N")

table2
########################################################
## Table 2 - Model 2 - First differences
########################################################

k <- 10000

# When government is litigant

m22.cond.govlit <- (data.m22[,"govislit"] == 1)
m22.data.govlit.X <- cbind(1,data.m22[m22.cond.govlit,vars.m22.X])

m22.govlit.X.set <- c()
m22.govlit.X.set <- my.setx(X = m22.data.govlit.X, X.set = m22.govlit.X.set, fun = my.mode, funs = list("percham" = mean, "CommObsPl" = modefunction, "CommObsDef" = modefunction, "CommIsPl" = modefunction, "CommIsDef" = modefunction))

m22.govlit.values <- my.values(list("normnetwobs" = c(0,0.1), "normwobsgov" = c(0,0.1)))
m22.govlit.Xs <- my.Xs(m22.govlit.X.set, m22.govlit.values)

m22.govlit.sim <- my.sim(m22.glm, k, m22.govlit.Xs)

m22.govlit.sim.fd <- m22.govlit.sim[,2] - m22.govlit.sim[,1]
mean(m22.govlit.sim.fd)
quantile(m22.govlit.sim.fd,0.025); quantile(m22.govlit.sim.fd,0.975)
mean(m22.govlit.sim[,1])
quantile(m22.govlit.sim[,1], 0.025); quantile(m22.govlit.sim[,1], 0.975)
mean(m22.govlit.sim[,2])
quantile(m22.govlit.sim[,2], 0.025); quantile(m22.govlit.sim[,2], 0.975)

# When government is not litigant

m22.cond.nogovlit <- (data.m22[,"govislit"] == 0)
m22.data.nogovlit.X <- cbind(1,data.m22[m22.cond.nogovlit,vars.m22.X])

# Set up representative data
m22.nogovlit.X.set <- c()
m22.nogovlit.X.set <- my.setx(X = m22.data.nogovlit.X, X.set = m22.nogovlit.X.set, fun = my.mode, funs = list("percham" = mean, "CommObsPl" = modefunction, "CommObsDef" = modefunction, "CommIsPl" = modefunction, "CommIsDef" = modefunction))
m22.nogovlit.values <- my.values(list("normnetwobs" = c(0,0.1)))

m22.nogovlit.Xs <- my.Xs(m22.nogovlit.X.set, m22.nogovlit.values)

# Simulation
m22.nogovlit.sim <- my.sim(m22.glm, k, m22.nogovlit.Xs)

# First Difference:
m22.nogovlit.sim.fd <- m22.nogovlit.sim[,2] - m22.nogovlit.sim[,1]
mean(m22.nogovlit.sim.fd)
quantile(m22.nogovlit.sim.fd, 0.025); quantile(m22.nogovlit.sim.fd, 0.975)
mean(m22.nogovlit.sim[,1])
quantile(m22.nogovlit.sim[,1], 0.025); quantile(m22.nogovlit.sim[,1], 0.975)
mean(m22.nogovlit.sim[,2])
quantile(m22.nogovlit.sim[,2], 0.025); quantile(m22.nogovlit.sim[,2], 0.975)



########################################################
## Table 3 - Model 1
########################################################

vars.m31.Y <- c("ECJPlAgree")
vars.m31.X <- c("normPl" , "normDef", "percham")
vars.m31.C <- c("article169")      # Filter condition
vars.m31 <- c(vars.m31.Y, vars.m31.X)
condition.m31 <- (data[,vars.m31.C] == 1)

data.m31 <- data[condition.m31,c(vars.m31, vars.G, vars.m31.C)]

formula.m31 <- as.formula(paste(c(vars.m31.Y,paste(unlist(vars.m31.X), collapse="+")), collapse="~"))

m31.glm <- glm(formula.m31, family=binomial(link = "probit"), data = data.m31, na.action = na.pass)
summary(m31.glm)
m31.my.se <- my.se(m31.glm, data.m31[,vars.G])
m31.my.se

infringe <- c(m31.my.se[2,1],m31.my.se[2,2],m31.my.se[3,1],m31.my.se[3,2],m31.my.se[4,1],m31.my.se[4,2],NA,NA,NA,NA,m31.my.se[1,1],m31.my.se[1,2],nrow(data.m31))
########################################################
## Table 3 - Model 2
########################################################

vars.m32.Y <- c("ECJPlAgree")
vars.m32.X <- c("normPl" , "normDef", "percham", "CommObsPl", "CommObsDef")
vars.m32.C <- c("article177", "govislit")      # Filter condition
vars.m32 <- c(vars.m32.Y, vars.m32.X)
condition.m32 <- ((data[,vars.m32.C[1]] == 1) & (data[,vars.m32.C[2]] == 0))

data.m32 <- data[condition.m32,c(vars.m32, vars.G, vars.m32.C)]
data.m32 <- na.omit(data.m32)

formula.m32 <- as.formula(paste(c(vars.m32.Y,paste(unlist(vars.m32.X), collapse="+")), collapse="~"))

m32.glm <- glm(formula.m32, family=binomial(link = "probit"), data = data.m32, na.action = na.pass)
summary(m32.glm)
m32.my.se <- my.se2(m32.glm, data.m32[,vars.G])
m32.my.se


art234 <- c(m32.my.se[2,1],m32.my.se[2,2],m32.my.se[3,1],m32.my.se[3,2],m32.my.se[4,1],m32.my.se[4,2],m32.my.se[5,1],m32.my.se[5,2],m32.my.se[6,1],m32.my.se[6,2],m32.my.se[1,1],m32.my.se[1,2], nrow(data.m32))
table3 <- cbind(infringe, art234)
rownames(table3) <- c("Total weighted observations for plaintiff", "se1","Total weighted observations for defendant","se2", "Chamber size", "se3", "Commission observation for plaintiff", "se4", "Commission observation for defendant", "se5","Constant", "se6","N")
colnames(table3) <- c("Infringement Proceeding", "Article 234 Case Without a Government Litigant")
table3

########################################################
## Table 3 - Model 1 - First differences
########################################################

k <- 10000
m31.cond.art169 <- (data.m31[,"article169"] == 1)
m31.data.art169.X <- cbind(1,data.m31[m31.cond.art169,vars.m31.X])

m31.art169.X.set <- c()
m31.art169.X.set <- my.setx(X = m31.data.art169.X, X.set = m31.art169.X.set, fun = my.mode, funs = list("percham" = median))
m31.art169.values <- my.values(list("normPl" = c(0,0.1))
m31.art169.Xs <- my.Xs(m31.art169.X.set, m31.art169.values)

m31.art169.sim <- my.sim(m31.glm, k, m31.art169.Xs)

m31.art169.sim.fd <- m31.art169.sim[,2] - m31.art169.sim[,1]
mean(m31.art169.sim.fd)
quantile(m31.art169.sim.fd,0.025);quantile(m31.art169.sim.fd,0.975)
mean(m31.art169.sim[,1])
quantile(m31.art169.sim[,1], 0.025); quantile(m31.art169.sim[,1], 0.975)
mean(m31.art169.sim[,2])
quantile(m31.art169.sim[,2], 0.025); quantile(m31.art169.sim[,2], 0.975)

########################################################
## Table 3 - Model 2 - First Differences
########################################################

k <- 10000
m32.cond.noart169 <- (data.m32[,"article177"] == 1) & (data.m32[,"govislit"] == 0)
m32.data.noart169.X <- cbind(1,data.m32[m32.cond.noart169,vars.m32.X])

m32.noart169.X.set <- c()
m32.noart169.X.set <- my.setx(X = m32.data.noart169.X, X.set = m32.noart169.X.set, fun = my.mode, funs = list("percham" = mean, "CommObsPl" = modefunction, "CommObsDef" = modefunction))
m32.noart169.values <- my.values(list("normPl" = c(0,0.1)))
m32.noart169.Xs <- my.Xs(m32.noart169.X.set, m32.noart169.values)

m32.noart169.sim <- my.sim(m32.glm, k, m32.noart169.Xs)

m32.noart169.sim.fd <- m32.noart169.sim[,2] - m32.noart169.sim[,1]
mean(m32.noart169.sim.fd)
quantile(m32.noart169.sim.fd, 0.025); quantile(m32.noart169.sim.fd, 0.975) 
mean(m32.noart169.sim[,1])
quantile(m32.noart169.sim[,1], 0.025); quantile(m32.noart169.sim[,1], 0.975)
mean(m32.noart169.sim[,2])
quantile(m32.noart169.sim[,2], 0.025); quantile(m32.noart169.sim[,2], 0.975)


#What about the difference in first differences?
diff.in.fd <- m31.art169.sim.fd - m32.noart169.sim.fd
mean(diff.in.fd)
quantile(diff.in.fd,0.025);quantile(diff.in.fd, 0.975)
########################################################
## Table 2 - Model 2 - Figure 1
########################################################

k <- 100

# Government is litigant
m22.f1.cond.govlit <- (data.m22[,"govislit"] == 1)
m22.f1.data.govlit.X <- cbind(1,data.m22[m22.f1.cond.govlit,vars.m22.X])

m22.f1.govlit.X.set <- c()
m22.f1.govlit.X.set <- my.setx(X = m22.f1.data.govlit.X, X.set = m22.f1.govlit.X.set, fun = my.mode, funs = list("percham" = median, "CommObsPl" = modefunction, "CommObsDef" = modefunction))

m22.f1.govlit.rng <- seq(from=-0.5, to = 0.5, by = 0.01)
m22.f1.govlit.values <- my.values(list("normnetwobs" = c(m22.f1.govlit.rng), "normwobsgov" = c(m22.f1.govlit.rng)))

m22.f1.govlit.Xs <- my.Xs(m22.f1.govlit.X.set, m22.f1.govlit.values)
m22.f1.govlit.sim <- my.sim(m22.glm, k, m22.f1.govlit.Xs)

# Government is not litigant
m22.f1.cond.nogovlit <- (data.m22[,"govislit"] == 0)
m22.f1.data.nogovlit.X <- cbind(1,data.m22[m22.f1.cond.nogovlit,vars.m22.X])

m22.f1.nogovlit.X.set <- c()
m22.f1.nogovlit.X.set <- my.setx(X = m22.f1.data.nogovlit.X, X.set = m22.f1.nogovlit.X.set, fun = my.mode, funs = list("percham" = median, "CommObsPl" = modefunction, "CommObsDef" = modefunction))

m22.f1.nogovlit.rng <- seq(from=-0.5, to = 0.5, by = 0.01)
m22.f1.nogovlit.values <- my.values(list("normnetwobs" = c(m22.f1.nogovlit.rng)))

m22.f1.nogovlit.Xs <- my.Xs(m22.f1.nogovlit.X.set, m22.f1.nogovlit.values)
m22.f1.nogovlit.sim <- my.sim(m22.glm, k, m22.f1.nogovlit.Xs)

## Plot where government is not litigant

plot(m22.f1.nogovlit.rng, colMeans(m22.f1.nogovlit.sim), ylim=c(0,1), col="black", pch=20, ylab="Probability of Ruling for Plaintiff", xlab="Net number of weighted government observations \n (as portion of all governments) for the plaintiff", xlim=c(m22.f1.nogovlit.rng[1],m22.f1.nogovlit.rng[length(m22.f1.nogovlit.rng)]), main = "The Impact of Government Observations on ECJ Ruling")
points(m22.f1.nogovlit.rng, apply(m22.f1.nogovlit.sim,2,quantile, 0.025), pch=20, cex=.2, col="black")
points(m22.f1.nogovlit.rng, apply(m22.f1.nogovlit.sim,2,quantile, 0.975), pch=20, cex=.2, col="black")

## Plot where government is litigant

points(m22.f1.govlit.rng, colMeans(m22.f1.govlit.sim), pch=20, col="gray")
points(m22.f1.govlit.rng, apply(m22.f1.govlit.sim,2,quantile, 0.025), cex=.2, pch=1, col="gray")
points(m22.f1.govlit.rng, apply(m22.f1.govlit.sim,2,quantile, 0.975), cex=.2, pch=1, col="gray")

legend(0,0.2, c("Government not litigant", "Government litigant"), col = c("black","gray"),
       text.col = 1, lty = c(1, 1), pch = c(20, 20),
       merge = TRUE)



########################################################
## Table 4 - Model 1
########################################################

vars.m41.Y <- c("ECJPlAgree")
vars.m41.X <- c("normnetwobs" , "govislit", "normwobsgov", "article177", "govlit177", "normwobs177", "normwobsgov177", "percham" , "CommIsPl" , "CommIsDef" , "CommObsPl", "CommObsDef")
vars.m41 <- c(vars.m41.Y, vars.m41.X)

data.m41 <- data[,c(vars.m41, vars.G)]
data.m41 <- na.omit(data.m41)

formula.m41 <- as.formula(paste(c(vars.m41.Y,paste(unlist(vars.m41.X), collapse="+")), collapse="~"))

m41.glm <- glm(formula.m41, family=binomial(link = "probit"), data = data.m41, na.action = na.pass)
summary(m41.glm)
m41.my.se <- my.se(m41.glm, data.m41[,vars.G])
m41.my.se


t4m1 <- c(m41.my.se[2,1],m41.my.se[2,2],m41.my.se[4,1],m41.my.se[4,2],m41.my.se[3,1],m41.my.se[3,2],m41.my.se[5,1],m41.my.se[5,2],m41.my.se[6,1],m41.my.se[6,2],m41.my.se[7,1],m41.my.se[7,2],m41.my.se[8,1],m41.my.se[8,2],m41.my.se[9,1],m41.my.se[9,2],m41.my.se[10,1],m41.my.se[10,2],m41.my.se[11,1],m41.my.se[11,2],m41.my.se[12,1],m41.my.se[12,2],m41.my.se[13,1],m41.my.se[13,2],NA,NA,m41.my.se[1,1],m41.my.se[1,2],nrow(data.m41))

########################################################
## Table 4 - Model 2
########################################################

vars.m42.Y <- c("ECJPlAgree")
vars.m42.X <- c("normnetwobs" , "govislit", "normwobsgov", "article177", "govlit177", "normwobs177", "normwobsgov177", "percham" , "CommIsPl" , "CommIsDef" , "CommObsPl", "CommObsDef", "AGforPl")
vars.m42 <- c(vars.m42.Y, vars.m42.X)

data.m42 <- data[,c(vars.m42, vars.G)]
data.m42 <- na.omit(data.m42)

formula.m42 <- as.formula(paste(c(vars.m42.Y,paste(unlist(vars.m42.X), collapse="+")), collapse="~"))

m42.glm <- glm(formula.m42, family=binomial(link = "probit"), data = data.m42, na.action = na.pass)
summary(m42.glm)
m42.my.se <- my.se(m42.glm, data.m42[,vars.G])
m42.my.se


t4m2 <- c(m42.my.se[2,1], m42.my.se[2,2],m42.my.se[4,1], m42.my.se[4,2],m42.my.se[3,1], m42.my.se[3,2],m42.my.se[5,1], m42.my.se[5,2],m42.my.se[6,1], m42.my.se[6,2],m42.my.se[7,1], m42.my.se[7,2],m42.my.se[8,1], m42.my.se[8,2],m42.my.se[9,1], m42.my.se[9,2],m42.my.se[10,1], m42.my.se[10,2],m42.my.se[11,1], m42.my.se[11,2],m42.my.se[12,1], m42.my.se[12,2],m42.my.se[13,1], m42.my.se[13,2],m42.my.se[14,1], m42.my.se[14,2],m42.my.se[1,1], m42.my.se[1,2],nrow(data.m42))

table4 <- cbind(t4m1, t4m2)
rownames(table4) <- c("Net weighted observations for plaintiff", "se1", "Net weighted observations for plaintiff x Government litigant","se2","Government litigant","se3","Article 234 case", "se4","Article 234 case x Government litigant","se5", "Net weighted observations for plaintiff x Article 234 case","se6","Net weighted observations for plaintiff x Article 234 case x Government litigant","se7", "Chamber size","se8", "Commission is plaintiff","se9", "Commission is defendant","se10","Commission observation for plaintiff","se11", "Commission observation for defendant","se12", "Advocate General for plaintiff","se13", "Constant","se14","N")
colnames(table4) <- c("Model 1", "Model 2")
table4
########################################################
## Table 4 - Model 1 - Figure 2
########################################################

k <- 100

# Preliminary Reference
m41.f2.cond.art177 <- (data.m41[,"article177"] == 1) & (data.m41[,"govislit"] == 1)
m41.f2.data.art177.X <- cbind(1,data.m41[m41.f2.cond.art177,vars.m41.X])

m41.f2.art177.X.set <- c()
m41.f2.art177.X.set <- my.setx(X = m41.f2.data.art177.X, X.set = m41.f2.art177.X.set, fun = my.mode, funs = list("percham" = median, "CommObsPl" = mean, "CommObsDef" = mean))

m41.f2.art177.rng <- seq(from=-0.5, to = 0.5, by = 0.01)
m41.f2.art177.values <- my.values(list("normnetwobs" = c(m41.f2.art177.rng), "normwobsgov" = c(m41.f2.art177.rng), "normwobs177" = c(m41.f2.art177.rng), "normwobsgov177" = c(m41.f2.art177.rng)))

m41.f2.art177.Xs <- my.Xs(m41.f2.art177.X.set, m41.f2.art177.values)
m41.f2.art177.sim <- my.sim(m41.glm, k, m41.f2.art177.Xs)


# No Preliminary Reference
m41.f2.cond.noart177 <- (data.m41[,"article177"] == 0) & (data.m41[,"govislit"] == 1)
m41.f2.data.noart177.X <- cbind(1,data.m41[m41.f2.cond.noart177,vars.m41.X])

m41.f2.noart177.X.set <- c()
m41.f2.noart177.X.set <- my.setx(X = m41.f2.data.noart177.X, X.set = m41.f2.noart177.X.set, fun = my.mode, funs = list("percham" = median, "CommObsPl" = mean, "CommObsDef" = mean))

m41.f2.noart177.rng <- seq(from=-0.5, to = 0.5, by = 0.01)
m41.f2.noart177.values <- my.values(list("normnetwobs" = c(m41.f2.noart177.rng), "normwobsgov" = c(m41.f2.noart177.rng)))

m41.f2.noart177.Xs <- my.Xs(m41.f2.noart177.X.set, m41.f2.noart177.values)
m41.f2.noart177.sim <- my.sim(m41.glm, k, m41.f2.noart177.Xs)


## Plot with non-preliminary reference cases

plot(m41.f2.noart177.rng, colMeans(m41.f2.noart177.sim), ylim=c(0,1), col="gray", pch=20, ylab="", xlab="Net number of weighted government observations \n (as portion of all governments) for the plaintiff", main = "Impact of Government Observations, Government as Litigant",xlim=c(m41.f2.noart177.rng[1],m41.f2.noart177.rng[length(m41.f2.noart177.rng)]))
points(m41.f2.noart177.rng, apply(m41.f2.noart177.sim,2,quantile, 0.025), pch=20, cex=.35, col="gray")
points(m41.f2.noart177.rng, apply(m41.f2.noart177.sim,2,quantile, 0.975), pch=20, cex=.35, col="gray")

## Plot with preliminary reference cases

points(m41.f2.art177.rng, colMeans(m41.f2.art177.sim), pch=20, col="black")
points(m41.f2.art177.rng, apply(m41.f2.art177.sim,2,quantile, 0.025), cex=.35, pch=1, col="black")
points(m41.f2.art177.rng, apply(m41.f2.art177.sim,2,quantile, 0.975), cex=.35, pch=1, col="black")

legend(0,0.2, c("Not Preliminary Reference", "Preliminary Reference"), col = c("gray","black"),
       text.col = 1, lty = c(1, 1), pch = c(20, 20),
       merge = TRUE)



########################################################
## Table 5 - Model 1 - Quantities of Interest
########################################################

## Main covariate
vars.m41.M <- c("normnetwobs")

## Relevant double interaction
vars.m41.D1 <- c("normwobsgov")
vars.m41.D2 <- c("normwobs177")

## Triple interaction
vars.m41.T <- c("normwobsgov177")

## Interaction vars
vars.m41.11 <- c(vars.m41.M, vars.m41.D1, vars.m41.D2, vars.m41.T)
vars.m41.10 <- c(vars.m41.M, vars.m41.D1)
vars.m41.01 <- c(vars.m41.M, vars.m41.D2)
vars.m41.00 <- c(vars.m41.M)

## Quantities of interest
zeta.m41.11 <- m41.glm$coefficients[vars.m41.M] + m41.glm$coefficients[vars.m41.D1] + m41.glm$coefficients[vars.m41.D2] + m41.glm$coefficients[vars.m41.T]
zeta.m41.10 <- m41.glm$coefficients[vars.m41.M] + m41.glm$coefficients[vars.m41.D1]
zeta.m41.01 <- m41.glm$coefficients[vars.m41.M] + m41.glm$coefficients[vars.m41.D2]
zeta.m41.00 <- m41.glm$coefficients[vars.m41.M]


## Standard errors

## Gradient

## Quantities of interest (see notes for more details):
## eta_11: b1 + b12 + b13 + b123
## eta_10: b1 + b12
## eta_01: b1 + b13
## eta_00: b1

## Order in which betas appear in the vcov matrix:
## b0 b1 b2 b12 b3 b23 b13 b123 [control correlates]

## Gradients
##  Grad11: <0  1  0  1   0  0   1   1>
##  Grad10: <0  1  0  1   0  0   0   0>
##  Grad01: <0  1  0  0   0  0   1   0>
##  Grad00: <1  0  0  0   0  0   0   0>   

## Retrieve HC variance-covariance matrix from probit fit
vcov.m41 <- my.vcov(m41.glm, data.m41[,vars.G])

grad.m41.11.eta <- c(0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0)
grad.m41.10.eta <- c(0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0)
grad.m41.01.eta <- c(0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0)
grad.m41.00.eta <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

se.m41.11 = sqrt(grad.m41.11.eta %*% vcov.m41 %*% grad.m41.11.eta)
se.m41.10 = sqrt(grad.m41.10.eta %*% vcov.m41 %*% grad.m41.10.eta)
se.m41.01 = sqrt(grad.m41.01.eta %*% vcov.m41 %*% grad.m41.01.eta)
se.m41.00 = sqrt(grad.m41.00.eta %*% vcov.m41 %*% grad.m41.00.eta)

## Put in a vector for convenience
t5m1 <- c(zeta.m41.11,se.m41.11,zeta.m41.10,se.m41.10,zeta.m41.01,se.m41.01,zeta.m41.00,se.m41.00)


########################################################
## Table 5 - Model 2 - Quantities of Interest
########################################################

vars.m42.M <- c("normnetwobs")
vars.m42.D1 <- c("normwobsgov")
vars.m42.D2 <- c("normwobs177")
vars.m42.T <- c("normwobsgov177")

zeta.m42.11 <- m42.glm$coefficients[vars.m42.M] + m42.glm$coefficients[vars.m42.D1] + m42.glm$coefficients[vars.m42.D2] + m42.glm$coefficients[vars.m42.T]
zeta.m42.10 <- m42.glm$coefficients[vars.m42.M] + m42.glm$coefficients[vars.m42.D1]
zeta.m42.01 <- m42.glm$coefficients[vars.m42.M] + m42.glm$coefficients[vars.m42.D2]
zeta.m42.00 <- m42.glm$coefficients[vars.m42.M]


vcov.m42 <- my.vcov(m42.glm, data.m42[,vars.G])

grad.m42.11.eta <- c(0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0)
grad.m42.10.eta <- c(0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
grad.m42.01.eta <- c(0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0)
grad.m42.00.eta <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

se.m42.11 = sqrt(grad.m42.11.eta %*% vcov.m42 %*% grad.m42.11.eta)
se.m42.10 = sqrt(grad.m42.10.eta %*% vcov.m42 %*% grad.m42.10.eta)
se.m42.01 = sqrt(grad.m42.01.eta %*% vcov.m42 %*% grad.m42.01.eta)
se.m42.00 = sqrt(grad.m42.00.eta %*% vcov.m42 %*% grad.m42.00.eta)


t5m2 <- c(zeta.m42.11,se.m42.11,zeta.m42.10,se.m42.10,zeta.m42.01,se.m42.01,zeta.m42.00,se.m42.00)
table5 <- cbind(t5m1, t5m2)
rownames(table5) <- c("Government is litigant and 234 case","se1","Government is litigant and not 234 case", "se2","Government is not litigant and 234 case", "se3","Government is not litigant and not 234 case","se4")
colnames(table5) <- c("Model 1", "Model 2")
table5
########################################################
## Table 6 - Model 1
########################################################

vars.m61.Y <- c("ECJPlAgree")
vars.m61.X <- c("perobs" , "govislit", "perobsgovlit", "percham" , "CommIsPl" , "CommIsDef" , "CommObsPl", "CommObsDef")
vars.m61 <- c(vars.m61.Y, vars.m61.X)

data.m61 <- data[,c(vars.m61, vars.G)]

formula.m61 <- as.formula(paste(c(vars.m61.Y,paste(unlist(vars.m61.X), collapse="+")), collapse="~"))

m61.glm <- glm(formula.m61, family=binomial(link = "probit"), data = data.m61, na.action = na.pass)
summary(m61.glm)
m61.my.se <- my.se(m61.glm, data.m61[,vars.G])
m61.my.se


########################################################
## Table 6 - Model 2
########################################################

vars.m62.Y <- c("ECJPlAgree")
vars.m62.X <- c("netGDPwobs" , "govislit", "netGDPwobsgov", "percham" , "CommIsPl" , "CommIsDef" , "CommObsPl", "CommObsDef")
vars.m62 <- c(vars.m62.Y, vars.m62.X)

data.m62 <- data[,c(vars.m62, vars.G)]

formula.m62 <- as.formula(paste(c(vars.m62.Y,paste(unlist(vars.m62.X), collapse="+")), collapse="~"))

m62.glm <- glm(formula.m62, family=binomial(link = "probit"), data = data.m62, na.action = na.pass)
summary(m62.glm)
m62.my.se <- my.se(m62.glm, data.m62[,vars.G])
m62.my.se


########################################################
## Table 6 - Model 3
########################################################

vars.m63.Y <- c("ECJPlAgree")
vars.m63.X <- c("nettradewobs" , "govislit", "nettradewobsgov", "percham" , "CommIsPl" , "CommIsDef" , "CommObsPl", "CommObsDef")
vars.m63 <- c(vars.m63.Y, vars.m63.X)

data.m63 <- data[,c(vars.m63, vars.G)]

formula.m63 <- as.formula(paste(c(vars.m63.Y,paste(unlist(vars.m63.X), collapse="+")), collapse="~"))

m63.glm <- glm(formula.m63, family=binomial(link = "probit"), data = data.m63, na.action = na.pass)
summary(m63.glm)
m63.my.se <- my.se(m63.glm, data.m63[,vars.G])
m63.my.se


########################################################
## Table 7 - Model 1
########################################################

vars.m71.Y <- c("ECJPlAgree")
vars.m71.X <- c("perobsPl" , "perobsDef", "percham")
vars.m71.C <- c("article169")      # Filter condition
vars.m71 <- c(vars.m71.Y, vars.m71.X)
condition.m71 <- (data[,vars.m71.C] == 1)

data.m71 <- data[condition.m71,c(vars.m71, vars.G, vars.m71.C)]

formula.m71 <- as.formula(paste(c(vars.m71.Y,paste(unlist(vars.m71.X), collapse="+")), collapse="~"))

m71.glm <- glm(formula.m71, family=binomial(link = "probit"), data = data.m71, na.action = na.pass)
summary(m71.glm)
m71.my.se <- my.se(m71.glm, data.m71[,vars.G])
m71.my.se


########################################################
## Table 7 - Model 2
########################################################

vars.m72.Y <- c("ECJPlAgree")
vars.m72.X <- c("GDPwobsPl" , "GDPwobsDef", "percham")
vars.m72.C <- c("article169")      # Filter condition
vars.m72 <- c(vars.m72.Y, vars.m72.X)
condition.m72 <- (data[,vars.m72.C] == 1)

data.m72 <- data[condition.m72,c(vars.m72, vars.G, vars.m72.C)]

formula.m72 <- as.formula(paste(c(vars.m72.Y,paste(unlist(vars.m72.X), collapse="+")), collapse="~"))

m72.glm <- glm(formula.m72, family=binomial(link = "probit"), data = data.m72, na.action = na.pass)
summary(m72.glm)
m72.my.se <- my.se(m72.glm, data.m72[,vars.G])
m72.my.se


########################################################
## Table 7 - Model 3
########################################################

vars.m73.Y <- c("ECJPlAgree")
vars.m73.X <- c("tradewobsPl" , "tradewobsDef", "percham")
vars.m73.C <- c("article169")      # Filter condition
vars.m73 <- c(vars.m73.Y, vars.m73.X)
condition.m73 <- (data[,vars.m73.C] == 1)

data.m73 <- data[condition.m73,c(vars.m73, vars.G, vars.m73.C)]

formula.m73 <- as.formula(paste(c(vars.m73.Y,paste(unlist(vars.m73.X), collapse="+")), collapse="~"))

m73.glm <- glm(formula.m73, family=binomial(link = "probit"), data = data.m73, na.action = na.pass)
summary(m73.glm)
m73.my.se <- my.se(m73.glm, data.m73[,vars.G])
m73.my.se


########################################################
## Table 8 - Model 1
########################################################

vars.m81.Y <- c("ECJPlAgree")
vars.m81.X <- c("perobsPl" , "perobsDef", "percham", "CommObsPl", "CommObsDef")
vars.m81.C <- c("article177", "govislit")      # Filter condition
vars.m81 <- c(vars.m81.Y, vars.m81.X)
condition.m81 <- ((data[,vars.m81.C[1]] == 1) & (data[,vars.m81.C[2]] == 0))

data.m81 <- data[condition.m81,c(vars.m81, vars.G, vars.m81.C)]
data.m81 <- na.omit(data.m81)

formula.m81 <- as.formula(paste(c(vars.m81.Y,paste(unlist(vars.m81.X), collapse="+")), collapse="~"))

m81.glm <- glm(formula.m81, family=binomial(link = "probit"), data = data.m81, na.action = na.pass)
summary(m81.glm)
m81.my.se <- my.se(m81.glm, data.m81[,vars.G])


########################################################
## Table 8 - Model 2
########################################################

vars.m82.Y <- c("ECJPlAgree")
vars.m82.X <- c("GDPwobsPl" , "GDPwobsDef", "percham", "CommObsPl", "CommObsDef")
vars.m82.C <- c("article177", "govislit")      # Filter condition
vars.m82 <- c(vars.m82.Y, vars.m82.X)
condition.m82 <- ((data[,vars.m82.C[1]] == 1) & (data[,vars.m82.C[2]] == 0))

data.m82 <- data[condition.m82,c(vars.m82, vars.G, vars.m82.C)]
data.m82 <- na.omit(data.m82)

formula.m82 <- as.formula(paste(c(vars.m82.Y,paste(unlist(vars.m82.X), collapse="+")), collapse="~"))

m82.glm <- glm(formula.m82, family=binomial(link = "probit"), data = data.m82, na.action = na.pass)
summary(m82.glm)
m82.my.se <- my.se(m82.glm, data.m82[,vars.G])


########################################################
## Table 8 - Model 3
########################################################

vars.m83.Y <- c("ECJPlAgree")
vars.m83.X <- c("tradewobsPl" , "tradewobsDef", "percham", "CommObsPl", "CommObsDef")
vars.m83.C <- c("article177", "govislit")      # Filter condition
vars.m83 <- c(vars.m83.Y, vars.m83.X)
condition.m83 <- ((data[,vars.m83.C[1]] == 1) & (data[,vars.m83.C[2]] == 0))

data.m83 <- data[condition.m83,c(vars.m83, vars.G, vars.m83.C)]
data.m83 <- na.omit(data.m83)

formula.m83 <- as.formula(paste(c(vars.m83.Y,paste(unlist(vars.m83.X), collapse="+")), collapse="~"))

m83.glm <- glm(formula.m83, family=binomial(link = "probit"), data = data.m83, na.action = na.pass)
summary(m83.glm)
m83.my.se <- my.se(m83.glm, data.m83[,vars.G])

























































data$CommObsDef[1197] <- 0

vars.G <- c("casenumber")
## Extract a robust, clustered variance-covariance matrix         
my.vcov <- function(fm, cluster){
         M <- length(unique(cluster))   
         N <- length(cluster)           
         K <- fm$rank                        
         dfc <- (M/(M-1))*((N-1)/(N-K))
         W <- as.vector(fm$prior.weights)
         # dfc <- M / (M - 1) // Another possible finite sample adjustment    
         uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
         vcovCL <- dfc*sandwich(fm, meat = crossprod(uj)/N)
         }
         

## Implements clustered vcov by default
my.vcov.c <- function(fm) my.vcov(fm, fm$data[,vars.G])

## Implements robust standard errors (people.su.se/~ma/clmclx.R)
my.se <- function(fm, cluster) coeftest(fm, my.vcov(fm, cluster))

## Implements robust standard errors (people.su.se/~ma/clmclx.R)
my.se.c <- function(fm) coeftest(fm, my.vcov.c(fm))


library(cem)
library(MatchIt)







#Table 2, Model 2
#let's just run the model real quick
vars.m22.Y <- c("ECJPlAgree")
vars.m22.X <- c("normnetwobs" , "govislit", "normwobsgov", "percham" , "CommIsPl" , "CommIsDef" , "CommObsPl", "CommObsDef")
vars.m22 <- c(vars.m22.Y, vars.m22.X)

data.m22 <- data[,c(vars.m22, vars.G)]

formula.m22 <- as.formula(paste(c(vars.m22.Y,paste(unlist(vars.m22.X), collapse="+")), collapse="~"))

m22.glm <- glm(formula.m22, family=binomial(link = "probit"), data = data.m22, na.action = na.pass)
summary(m22.glm)
m22.my.se <- my.se(m22.glm, data.m22[,vars.G])
m22.my.se

#calculating imbalance
imbal1 <- imbalance(data$normnetwobs,data = data.m22,drop = c("casenumber"))
imbal1


#let's do the convex hull analysis
library(WhatIf)

#for the litigants: 50 out of 101 fall within the hull:-0.22 to 0.28
cfactmatlitig <- matrix(data = NA, nrow = 101, ncol = 8)
k <- -0.5
for(i in 1:101){
	cfactmatlitig[i,] <- c(k, 1, k,median(data$percham), modefunction(data$CommIsPl),modefunction(data$CommIsDef), modefunction(data$CommObsPl), modefunction(data$CommObsDef))
	k <- k + 0.01
	}	
chullq <- whatif(data = m22.glm, cfact = cfactmatlitig)
chullq$in.hull
hullresultlitig <- cbind(chullq$in.hull, cfactmatlitig[,1])


#for non-litigants: 64 of 101 fall within the convex hull: -0.44 to 0.19
cfactmatnonlitig <- matrix(data = NA, nrow = 101, ncol = 8)
k <- -0.5
for(i in 1:101){
	cfactmatnonlitig[i,] <- c(k, 0, 0,median(data$percham), modefunction(data$CommIsPl),modefunction(data$CommIsDef), modefunction(data$CommObsPl), modefunction(data$CommObsDef))
	k <- k + 0.01
	}	
chullq2 <- whatif(data = m22.glm, cfact = cfactmatnonlitig)
chullq2$in.hull
hullresultnonlitig <- cbind(chullq2$in.hull, cfactmatnonlitig[,1])


#let's generate the figure
k <- 100

# Government is litigant
m22.f1.cond.govlit <- (data.m22[,"govislit"] == 1)
m22.f1.data.govlit.X <- cbind(1,data.m22[m22.f1.cond.govlit,vars.m22.X])

m22.f1.govlit.X.set <- c()
m22.f1.govlit.X.set <- my.setx(X = m22.f1.data.govlit.X, X.set = m22.f1.govlit.X.set, fun = my.mode, funs = list("percham" = median, "CommObsPl" = modefunction, "CommObsDef" = modefunction))

m22.f1.govlit.rng <- seq(from=-0.50, to = 0.50, by = 0.01)
m22.f1.govlit.values <- my.values(list("normnetwobs" = c(m22.f1.govlit.rng), "normwobsgov" = c(m22.f1.govlit.rng)))

m22.f1.govlit.Xs <- my.Xs(m22.f1.govlit.X.set, m22.f1.govlit.values)
m22.f1.govlit.sim <- my.sim(m22.glm, k, m22.f1.govlit.Xs)

# Government is not litigant
m22.f1.cond.nogovlit <- (data.m22[,"govislit"] == 0)
m22.f1.data.nogovlit.X <- cbind(1,data.m22[m22.f1.cond.nogovlit,vars.m22.X])

m22.f1.nogovlit.X.set <- c()
m22.f1.nogovlit.X.set <- my.setx(X = m22.f1.data.nogovlit.X, X.set = m22.f1.nogovlit.X.set, fun = my.mode, funs = list("percham" = median, "CommObsPl" = modefunction, "CommObsDef" = modefunction))

m22.f1.nogovlit.rng <- seq(from=-0.5, to = 0.5, by = 0.01)
m22.f1.nogovlit.values <- my.values(list("normnetwobs" = c(m22.f1.nogovlit.rng)))

m22.f1.nogovlit.Xs <- my.Xs(m22.f1.nogovlit.X.set, m22.f1.nogovlit.values)
m22.f1.nogovlit.sim <- my.sim(m22.glm, k, m22.f1.nogovlit.Xs)

## Plot where government is not litigant

plot(m22.f1.nogovlit.rng[1:6], colMeans(m22.f1.nogovlit.sim)[1:6], ylim=c(0,1), col="black", pch=24, ylab="Probability of Ruling for Plaintiff", xlab="Net number of weighted government observations \n (as portion of all governments) for the plaintiff", xlim=c(m22.f1.nogovlit.rng[1],m22.f1.nogovlit.rng[length(m22.f1.nogovlit.rng)]), main = "Figure 1, Carrubba et al (2008): The Impact \n of Government Observations on ECJ Ruling", cex = 1.5)
points(m22.f1.nogovlit.rng[7:70], colMeans(m22.f1.nogovlit.sim)[7:70], ylim=c(0,1), col="black", pch=20, cex = 1.5)
points(m22.f1.nogovlit.rng[71:101], colMeans(m22.f1.nogovlit.sim)[71:101], ylim=c(0,1), col="black", pch=24, cex = 1.5)
points(m22.f1.nogovlit.rng, apply(m22.f1.nogovlit.sim,2,quantile, 0.025), pch=20, cex=.2, col="black")
points(m22.f1.nogovlit.rng, apply(m22.f1.nogovlit.sim,2,quantile, 0.975), pch=20, cex=.2, col="black")

## Plot where government is litigant

points(m22.f1.govlit.rng[1:28], colMeans(m22.f1.govlit.sim)[1:28], pch=24, col="gray", cex = 1.5)
points(m22.f1.govlit.rng[29:78], colMeans(m22.f1.govlit.sim)[29:78], pch=20, col="gray", cex = 1.5)
points(m22.f1.govlit.rng[79:101], colMeans(m22.f1.govlit.sim)[79:101], pch=24, col="gray", cex = 1.5)
points(m22.f1.govlit.rng, apply(m22.f1.govlit.sim,2,quantile, 0.025), cex=.2, pch=1, col="gray")
points(m22.f1.govlit.rng, apply(m22.f1.govlit.sim,2,quantile, 0.975), cex=.2, pch=1, col="gray")

legend(0,0.2, c("Government not litigant", "Government litigant"), col = c("black","gray"),
       text.col = 1, lty = c(1,1),lwd = c(3,3), merge = TRUE)

legend(-0.4,0.9, c("Within convex hull","Outside convex hull"), pch = c(20,24), col = c("black","black"))












#take a quick look at the data
dataobs <- data[data$normnetwobs != 0,]
quantile(dataobs$normnetwobs, probs = seq(0,1,0.05))


#now let's do some coarsened exact matching to solve our pretty serious model dependence problem
data$treatmentcoarse <- ifelse(data$normnetwobs <= -0.50,0,ifelse(data$normnetwobs <= -0.25,1,ifelse(data$normnetwobs <= -0.13,2,ifelse(data$normnetwobs <= -0.09,3,ifelse(data$normnetwobs < 0,4,ifelse(data$normnetwobs == 0,5,ifelse(data$normnetwobs <= 0.09,6,ifelse(data$normnetwobs <= 0.18,7,8))))))))



mod2datacoarse <- cbind(data$ECJPlAgree,data$normnetwobs, data$normwobsgov, data$govislit, data$percham, data$CommIsPl, data$CommIsDef, data$CommObsPl, data$CommObsDef,data$treatmentcoarse,data$AGforPl)
mod2datacoarse <- data.frame(mod2datacoarse)
colnames(mod2datacoarse) <- c("ECJPlAgree","normnetwobs", "normwobsgov","govislit", "percham", "CommIsPl", "CommIsDef", "CommObsPl","CommObsDef", "treatmentcoarse","AGforPl")

perchamcut <- c(0.1,0.25,0.40,0.59,0.70,0.87)
user.match <- cem(treatment = "treatmentcoarse", data = mod2datacoarse, drop = c("ECJPlAgree","normwobsgov","normnetwobs"))
user.match
#weights matrix
weights <- user.match$w

#running our regression
vars.m22.Y <- c("ECJPlAgree")
vars.m22.X <- c("normnetwobs")
vars.m22 <- c(vars.m22.Y, vars.m22.X)

data.m22 <- data[,c(vars.m22, vars.G)]

formula.m22 <- as.formula(paste(c(vars.m22.Y,paste(unlist(vars.m22.X), collapse="+")), collapse="~"))

m22.glm <- glm(formula.m22, family=binomial(link = "probit"), data = data.m22, na.action = na.pass, weights = weights)
summary(m22.glm)
m22.my.se <- my.se(m22.glm, data.m22[,vars.G])
m22.my.se

#let's simulate
zout <- zelig(ECJPlAgree ~ normnetwobs, model = "probit", data = data.m22,weights = weights)
zsetx0 <- setx(zout, normnetwobs = 0)
zsetx1 <- setx(zout, normnetwobs = 0.1)
s.out0 <- sim(zout, x = zsetx0)
summary(s.out0)
s.out1 <- sim(zout, x = zsetx1)
summary(s.out1)
sdiff <- sim(zout, zsetx0, zsetx1)
summary(sdiff)



#test the convex hull!
cfactmatlitig <- matrix(data = NA, nrow = 101, ncol = 1)
k <- -0.5
for(i in 1:101){
	cfactmatlitig[i,] <- c(k)
	k <- k + 0.01
	}	
chullq <- whatif(data = m22.glm, cfact = cfactmatlitig)
chullq$in.hull
hullresultlitig <- cbind(chullq$in.hull, cfactmatlitig[,1])




#now, let's try to do this with govlit and normobs unmatched; these are bizarre results, particularly with the standard errors?
user.match <- cem(treatment = "treatmentcoarse", data = mod2datacoarse, drop = c("ECJPlAgree","normwobsgov","normnetwobs","govislit"))
user.match
#weights matrix
weights <- user.match$w

#running our regression
vars.m22.Y <- c("ECJPlAgree")
vars.m22.X <- c("normnetwobs","CommObsPl","govislit","normwobsgov")
vars.m22 <- c(vars.m22.Y, vars.m22.X)

data.m22 <- data[,c(vars.m22, vars.G)]

formula.m22 <- as.formula(paste(c(vars.m22.Y,paste(unlist(vars.m22.X), collapse="+")), collapse="~"))

m22.glm <- glm(formula.m22, family=binomial(link = "probit"), data = data.m22, na.action = na.pass, weights = weights)
summary(m22.glm)
m22.my.se <- my.se(m22.glm, data.m22[,vars.G])
m22.my.se

#let's do the simulation analysis
zout <- zelig(ECJPlAgree ~ normnetwobs + CommObsPl + govislit + normwobsgov, model = "probit",data = data.m22, weights = weights)
zsetx0 <- setx(zout, normnetwobs = 0, govislit = 1, normwobsgov = 0)
zsetx1 <- setx(zout, normnetwobs = 0.1, govislit = 1, normwobsgov = 0.1)
zsetx2 <- setx(zout, normnetwobs = 0, govislit = 0, normwobsgov = 0)
zsetx3 <- setx(zout, normnetwobs = 0.1, govislit = 0, normwobsgov = 0)
zsimx1 <- sim(zout, zsetx0, zsetx1)
summary(zsimx1)
zsim0 <- sim(zout, zsetx0)
summary(zsim0)
zsim1 <- sim(zout, zsetx1)
summary(zsim1)
zsimx3 <- sim(zout, zsetx2, zsetx3)
summary(zsimx3)

zout <- zelig(ECJPlAgree ~ normnetwobs, model = "probit", data = data.m22,weights = weights)
zsetx0 <- setx(zout, normnetwobs = 0)
zsetx1 <- setx(zout, normnetwobs = 0.1)
s.out0 <- sim(zout, x = zsetx0)
summary(s.out0)
s.out1 <- sim(zout, x = zsetx1)
summary(s.out1)
sdiff <- sim(zout, zsetx0, zsetx1)
summary(sdiff)



#test the convex hull
#test the convex hull!
cfactmatlitig <- matrix(data = NA, nrow = 101, ncol = 4)
k <- -0.5
for(i in 1:101){
	cfactmatlitig[i,] <- c(k,0,0,0)
	k <- k + 0.01
	}	
chullq <- whatif(data = m22.glm, cfact = cfactmatlitig)
chullq$in.hull
hullresultlitig <- cbind(chullq$in.hull, cfactmatlitig[,1])



#table 3: infringement proceedings
vars.m31.Y <- c("ECJPlAgree")
vars.m31.X <- c("normPl" , "normDef", "percham")
vars.m31.C <- c("article169")      # Filter condition
vars.m31 <- c(vars.m31.Y, vars.m31.X)
condition.m31 <- (data[,vars.m31.C] == 1)

data.m31 <- data[condition.m31,c(vars.m31, vars.G, vars.m31.C)]

formula.m31 <- as.formula(paste(c(vars.m31.Y,paste(unlist(vars.m31.X), collapse="+")), collapse="~"))

m31.glm <- glm(formula.m31, family=binomial(link = "probit"), data = data.m31, na.action = na.pass)
summary(m31.glm)
m31.my.se <- my.se(m31.glm, data.m31[,vars.G])
m31.my.se

imbal2 <- imbalance(data.m31$normPl,data = data.m31,drop = c("casenumber"))
imbal2

#now, let's do some convex hull analysis: if we hold constant defendant at 0, then only 0 to 0.13 observations in the plaintiff's favor are legit. This gets even worse if there's any defendant's observations (with 0.1 defendant observations, only up to 0.07 is legit)
cfactmatinfringe <- matrix(data = NA, nrow = 101, ncol = 3)
k <- 0
for(i in 1:101){
	cfactmatinfringe[i,] <- c(k,0,mean(data$percham))
	k <- k + 0.01
	}	
chullq <- whatif(data = m31.glm, cfact = cfactmatinfringe)
chullq$in.hull


cfactmatinfringe <- matrix(data = NA, nrow = 101, ncol = 3)
k <- 0
for(i in 1:101){
	cfactmatinfringe[i,] <- c(k,0.1,mean(data$percham))
	k <- k + 0.01
	}	
chullq <- whatif(data = m31.glm, cfact = cfactmatinfringe)
chullq$in.hull

#up to 0.17 defendant observations when plaintiff observations at 0; when plaintiff observations at 0.1, then only up to 0.06 -> that's some serious model dependence
cfactmatinfringe <- matrix(data = NA, nrow = 101, ncol = 3)
k <- 0
for(i in 1:101){
	cfactmatinfringe[i,] <- c(0,k,mean(data$percham))
	k <- k + 0.01
	}	
chullq <- whatif(data = m31.glm, cfact = cfactmatinfringe)
chullq$in.hull


cfactmatinfringe <- matrix(data = NA, nrow = 101, ncol = 3)
k <- 0
for(i in 1:101){
	cfactmatinfringe[i,] <- c(0.1,k,mean(data$percham))
	k <- k + 0.01
	}	
chullq <- whatif(data = m31.glm, cfact = cfactmatinfringe)
chullq$in.hull


#in how many cases were there any observations
datam31 <- data[data$article169 == 1,]
datam231 <- data[data$article169 == 0,]
datanormPlDef <- datam31[(datam31$normPl > 0) | (datam31$normDef > 0),]
datanormPl <- datam31[datam31$normPl > 0,]
datanormDef <- datam31[datam31$normDef > 0,]
datanormPlDef <- datam231[(datam231$normPl > 0) | (datam231$normDef > 0),]

treatnormPl <- ifelse(datam31$normPl == 0,0,ifelse(datam31$normPl >= 0.11,1,2))
mod3data <- cbind(datam31$ECJPlAgree,datam31$normPl, datam31$normDef, datam31$percham, datam31$CommObsPl, datam31$CommObsDef,datam31$AGforPl, treatnormPl, treatnormDef)
mod3data <- data.frame(mod3data)
colnames(mod3data) <- c("ECJPlAgree","normPl", "normDef","percham", "CommObsPl","CommObsDef","AGforPl","treatnormPl","treatnormDef")

perchamcut <- c(0.1,0.25,0.40,0.59,0.70,0.87)
perchamcut2 <- c(0.1,0.4,0.7)

user.match <- cem(treatment = "treatnormPl", data = mod3data,drop = c("ECJPlAgree","normPl", "CommObsPl","CommObsDef","treatnormDef","AGforPl"))
user.match
#weights matrix
weights <- user.match$w


#let's run the regression on this matched shit, for normPl - the effect is significant
vars.m31.Y <- c("ECJPlAgree")
vars.m31.X <- c("normPl","percham")
vars.m31.C <- c("article169")      # Filter condition
vars.m31 <- c(vars.m31.Y, vars.m31.X)
condition.m31 <- (data[,vars.m31.C] == 1)

data.m31 <- datam31[,c(vars.m31, vars.G)]

formula.m31 <- as.formula(paste(c(vars.m31.Y,paste(unlist(vars.m31.X), collapse="+")), collapse="~"))

m31.glm <- glm(formula.m31, family=binomial(link = "probit"), data = data.m31, na.action = na.pass, weights = weights)
summary(m31.glm)
m31.my.se <- my.se(m31.glm, data.m31[,vars.G])
m31.my.se

zout <- zelig(ECJPlAgree ~ normPl + percham, data = data.m31, model = "probit", weights = weights)
setx1 <- setx(zout, normPl = 0, percham = mean(data.m31$percham))
setx2 <- setx(zout, normPl = 0.1, percham = mean(data.m31$percham))
simdif <- sim(zout, setx1, setx2)

#Let's do observations for the defendant
treatnormDef <- ifelse(datam31$normDef == 0,0,ifelse(datam31$normDef <= 0.13,1,2))
mod3data <- cbind(datam31$ECJPlAgree,datam31$normPl, datam31$normDef, datam31$percham, datam31$CommObsPl, datam31$CommObsDef,datam31$AGforPl, treatnormPl, treatnormDef)
mod3data <- data.frame(mod3data)
colnames(mod3data) <- c("ECJPlAgree","normPl", "normDef","percham", "CommObsPl","CommObsDef","AGforPl","treatnormPl","treatnormDef")

perchamcut <- c(0.1,0.25,0.40,0.59,0.70,0.87)
perchamcut2 <- c(0.1,0.4,0.7)

user.match <- cem(treatment = "treatnormDef", data = mod3data, drop = c("ECJPlAgree","normDef","AGforPl","treatnormPl","CommObsPl","CommObsDef"))
user.match
#weights matrix
weights <- user.match$w


#let's run the regression on this matched shit, for normDef
vars.m31.Y <- c("ECJPlAgree")
vars.m31.X <- c("normPl","normDef","percham")
vars.m31.C <- c("article169")      # Filter condition
vars.m31 <- c(vars.m31.Y, vars.m31.X)
condition.m31 <- (data[,vars.m31.C] == 1)

data.m31 <- datam31[,c(vars.m31, vars.G)]

formula.m31 <- as.formula(paste(c(vars.m31.Y,paste(unlist(vars.m31.X), collapse="+")), collapse="~"))

m31.glm <- glm(formula.m31, family=binomial(link = "probit"), data = data.m31, na.action = na.pass, weights = weights)
summary(m31.glm)
m31.my.se <- my.se(m31.glm, data.m31[,vars.G])
m31.my.se

zout <- zelig(ECJPlAgree ~ normDef + normPl + percham, data = data.m31, model = "probit", weights = weights)
setx1 <- setx(zout, normDef = 0, normPl = 0, percham = mean(data.m31$percham))
setx2 <- setx(zout, normDef = 0.1, normPl = 0, percham = mean(data.m31$percham))
simdif <- sim(zout, setx1, setx2)
















#how about in article 234 cases without a government litigant
datamod32 <- data[(data$article177 == 1) & (data$govislit == 0),]
datamod32normPl <- datamod32[datamod32$normPl > 0,]
datamod32normDef <- datamod32[datamod32$normDef > 0,]

#let's run the original model
vars.m32.Y <- c("ECJPlAgree")
vars.m32.X <- c("normPl" , "normDef", "percham", "CommObsPl", "CommObsDef")
vars.m32.C <- c("article177", "govislit")      # Filter condition
vars.m32 <- c(vars.m32.Y, vars.m32.X)
condition.m32 <- ((data[,vars.m32.C[1]] == 1) & (data[,vars.m32.C[2]] == 0))

data.m32 <- data[condition.m32,c(vars.m32, vars.G, vars.m32.C)]
data.m32 <- na.omit(data.m32)

formula.m32 <- as.formula(paste(c(vars.m32.Y,paste(unlist(vars.m32.X), collapse="+")), collapse="~"))

m32.glm <- glm(formula.m32, family=binomial(link = "probit"), data = data.m32, na.action = na.pass)
summary(m32.glm)
m32.my.se <- my.se2(m32.glm, data.m32[,vars.G])
m32.my.se


imbal2 <- imbalance(data.m32$normPl,data = data.m31,drop = c("casenumber"))
imbal2



#for the plaintiff weighted observations, up to 0.51 if defendant is 0 (if Commission observation for Plaintiff); up to 0.21 if defendant is 0 (assuming no obs for Plaintiff) and 0.17 if defendant has 0.1
cfactmatinfringe <- matrix(data = NA, nrow = 101, ncol = 5)
k <- 0
for(i in 1:101){
	cfactmatinfringe[i,] <- c(k,0,mean(data$percham),0,0)
	k <- k + 0.01
	}	
chullq <- whatif(data = m32.glm, cfact = cfactmatinfringe)
chullq$in.hull


cfactmatinfringe <- matrix(data = NA, nrow = 101, ncol = 5)
k <- 0
for(i in 1:101){
	cfactmatinfringe[i,] <- c(k,0.1,mean(data$percham),0,0)
	k <- k + 0.01
	}	
chullq <- whatif(data = m32.glm, cfact = cfactmatinfringe)
chullq$in.hull

#up to 0.46 obs for defendant when all else at 0; up to 0.27 when defendant up to 0.1
cfactmatinfringe <- matrix(data = NA, nrow = 101, ncol = 5)
k <- 0
for(i in 1:101){
	cfactmatinfringe[i,] <- c(0,k,mean(data$percham),0,0)
	k <- k + 0.01
	}	
chullq <- whatif(data = m32.glm, cfact = cfactmatinfringe)
chullq$in.hull


cfactmatinfringe <- matrix(data = NA, nrow = 101, ncol = 5)
k <- 0
for(i in 1:101){
	cfactmatinfringe[i,] <- c(0.1,k,mean(data$percham),1,0)
	k <- k + 0.01
	}	
chullq <- whatif(data = m32.glm, cfact = cfactmatinfringe)
chullq$in.hull



#now, let's try some matching to see how this works; first for the plaintiff observations
treatnormPl32 <- ifelse(datamod32$normPl == 0, 0,ifelse(datamod32$normPl <= 0.12,1,ifelse(datamod32$normPl <= 0.15,2,ifelse(datamod32$normPl <= 0.3,3,4))))

mod32data <- cbind(datamod32$ECJPlAgree,datamod32$normPl, datamod32$normDef, datamod32$percham, datamod32$CommObsPl, datamod32$CommObsDef,datamod32$AGforPl, treatnormPl32)
mod32data <- data.frame(mod32data)
colnames(mod32data) <- c("ECJPlAgree","normPl", "normDef","percham", "CommObsPl","CommObsDef","AGforPl","treatnormPl32")

perchamcut <- c(0.1,0.25,0.40,0.59,0.70,0.87)
perchamcut2 <- c(0.1,0.4,0.7)

user.match <- cem(treatment = "treatnormPl32", data = mod32data, drop = c("ECJPlAgree","normPl","AGforPl"))
user.match
#weights matrix
weights <- user.match$w


#run the regressions
vars.m31.Y <- c("ECJPlAgree")
vars.m31.X <- c("normPl","normDef","percham","CommObsPl")
vars.m31.C <- c("article169")      # Filter condition
vars.m31 <- c(vars.m31.Y, vars.m31.X)
condition.m31 <- (data[,vars.m31.C] == 1)

data.m31 <- datamod32[,c(vars.m31, vars.G)]

formula.m31 <- as.formula(paste(c(vars.m31.Y,paste(unlist(vars.m31.X), collapse="+")), collapse="~"))

m31.glm <- glm(formula.m31, family=binomial(link = "probit"), data = data.m31, na.action = na.pass, weights = weights)
summary(m31.glm)
m31.my.se <- my.se(m31.glm, data.m31[,vars.G])
m31.my.se

zout <- zelig(ECJPlAgree ~ normDef + normPl + percham + CommObsPl, data = data.m32, model = "probit", weights = weights)
setx1 <- setx(zout, normPl = 0, normDef = 0, percham = mean(data.m32$percham), CommObsPl = 0)
setx2 <- setx(zout, normPl = 0.1, normDef = 0, percham = mean(data.m32$percham), CommObsPl = 0)
simdif <- sim(zout, setx1, setx2)






#how about the defendants observations?
treatnormDef32 <- ifelse(datamod32$normDef == 0, 0,ifelse(datamod32$normDef <= 0.1,1,ifelse(datamod32$normDef <= 0.14,2,ifelse(datamod32$normDef <= 0.3,3,4))))

mod32data <- cbind(datamod32$ECJPlAgree,datamod32$normPl, datamod32$normDef, datamod32$percham, datamod32$CommObsPl, datamod32$CommObsDef,datamod32$AGforPl, treatnormDef32)
mod32data <- data.frame(mod32data)
colnames(mod32data) <- c("ECJPlAgree","normPl", "normDef","percham", "CommObsPl","CommObsDef","AGforPl","treatnormDef32")

user.match <- cem(treatment = "treatnormDef32", data = mod32data, drop = c("ECJPlAgree","normDef","treatnorm","AGforPl"))
user.match
#weights matrix
weights <- user.match$w


#run the regressions
vars.m31.Y <- c("ECJPlAgree")
vars.m31.X <- c("normPl","normDef","percham","CommObsPl","CommObsDef")
vars.m31.C <- c("article169")      # Filter condition
vars.m31 <- c(vars.m31.Y, vars.m31.X)
condition.m31 <- (data[,vars.m31.C] == 1)

data.m31 <- datamod32[,c(vars.m31, vars.G)]

formula.m31 <- as.formula(paste(c(vars.m31.Y,paste(unlist(vars.m31.X), collapse="+")), collapse="~"))

m31.glm <- glm(formula.m31, family=binomial(link = "probit"), data = data.m31, na.action = na.pass, weights = weights)
summary(m31.glm)
m31.my.se <- my.se(m31.glm, data.m31[,vars.G])
m31.my.se

zout <- zelig(ECJPlAgree ~ normDef + normPl + percham + CommObsPl, data = data.m32, model = "probit", weights = weights)
setx1 <- setx(zout, normDef = 0, normPl = 0, percham = mean(data.m32$percham), CommObsPl = 0, CommObsDef = 0)
setx2 <- setx(zout, normDef = 0.1, normPl = 0, percham = mean(data.m32$percham), CommObsPl = 0, CommObsDef = 0)
simdif <- sim(zout, setx1, setx2)











#table 4: let's replicate it!
vars.m41.Y <- c("ECJPlAgree")
vars.m41.X <- c("normnetwobs" , "govislit", "normwobsgov", "article177", "govlit177", "normwobs177", "normwobsgov177", "percham" , "CommIsPl" , "CommIsDef" , "CommObsPl", "CommObsDef")
vars.m41 <- c(vars.m41.Y, vars.m41.X)

data.m41 <- data[,c(vars.m41, vars.G)]
data.m41 <- na.omit(data.m41)

formula.m41 <- as.formula(paste(c(vars.m41.Y,paste(unlist(vars.m41.X), collapse="+")), collapse="~"))

m41.glm <- glm(formula.m41, family=binomial(link = "probit"), data = data.m41, na.action = na.pass)
summary(m41.glm)
m41.my.se <- my.se(m41.glm, data.m41[,vars.G])
m41.my.se


#some convex hull analysis first; let's do Figure 2
#for article 234 government litigant, from -0.21 to 0.26
cfactmat41 <- matrix(data = NA, nrow = 101, ncol = 12)
k <- -0.5
for(i in 1:101){
	cfactmat41[i,] <- c(k,1,k,1,1,k,k,mean(data$percham),0,0,0,0)
	k <- k + 0.01
	}	
chullq <- whatif(data = m41.glm, cfact = cfactmat41)
chullq$in.hull

#for litigant non-234, none of the counterfactuals fall within the convex hull -> some serious extrapolation!
cfactmat41 <- matrix(data = NA, nrow = 101, ncol = 12)
k <- -0.5
for(i in 1:101){
	cfactmat41[i,] <- c(k,1,k,0,0,0,0,mean(data$percham),0,0,0,0)
	k <- k + 0.01
	}	
chullq <- whatif(data = m41.glm, cfact = cfactmat41)
chullq$in.hull

#for non-litigant 234, from -0.46 to 0.23
cfactmat41 <- matrix(data = NA, nrow = 101, ncol = 12)
k <- -0.5
for(i in 1:101){
	cfactmat41[i,] <- c(k,0,0,1,0,k,0,mean(data$percham),0,0,0,0)
	k <- k + 0.01
	}	
chullq <- whatif(data = m41.glm, cfact = cfactmat41)
chullq$in.hull

#for non-litigant non-234, only from -0.11 to 0 are we not extrapolating
cfactmat41 <- matrix(data = NA, nrow = 101, ncol = 12)
k <- -0.5
for(i in 1:101){
	cfactmat41[i,] <- c(k,0,0,0,0,0,0,mean(data$percham),0,0,0,0)
	k <- k + 0.01
	}	
chullq <- whatif(data = m41.glm, cfact = cfactmat41)
chullq$in.hull




#how about matching for the model including the AG?
data$treatmentcoarse <- ifelse(data$normnetwobs <= -0.25,0,ifelse(data$normnetwobs <= -0.13,1,ifelse(data$normnetwobs <= -0.07,2,ifelse(data$normnetwobs < 0,3,ifelse(data$normnetwobs == 0,4,ifelse(data$normnetwobs <= 0.10,5,6))))))

mod41datacoarse <- cbind(data$ECJPlAgree,data$normnetwobs,data$normwobsgov, data$govislit, data$article177, data$govlit177,data$normwobs177,data$normwobsgov177,data$percham, data$CommIsPl, data$CommIsDef, data$CommObsPl, data$CommObsDef,data$treatmentcoarse,data$AGforPl)
mod41datacoarse <- data.frame(mod41datacoarse)
colnames(mod41datacoarse) <- c("ECJPlAgree","normnetwobs", "normwobsgov","govislit", "article177","govlit177","normwobs177","normwobsgov177","percham", "CommIsPl", "CommIsDef", "CommObsPl","CommObsDef", "treatmentcoarse","AGforPl")


perchamcut <- c(0.1,0.25,0.40,0.59,0.70,0.87)
user.match <- cem(treatment = "treatmentcoarse", data = mod41datacoarse, drop = c("ECJPlAgree","normwobsgov","normnetwobs","normwobsgov","govislit","article177","govlit177","normwobs177","normwobsgov177","AGforPl","percham"))
user.match
#weights matrix
weights <- user.match$w

vars.m41.Y <- c("ECJPlAgree")
vars.m41.X <- c("normnetwobs" , "govislit", "normwobsgov", "article177", "govlit177", "normwobs177", "normwobsgov177", "percham" ,"CommObsPl","AGforPl")
vars.m41 <- c(vars.m41.Y, vars.m41.X)

data.m41 <- data[,c(vars.m41, vars.G)]
data.m41 <- na.omit(data.m41)

formula.m41 <- as.formula(paste(c(vars.m41.Y,paste(unlist(vars.m41.X), collapse="+")), collapse="~"))

m41.glm <- glm(formula.m41, family=binomial(link = "probit"), data = data.m41, na.action = na.pass, weights = weights)
summary(m41.glm)
m41.my.se <- my.se(m41.glm, data.m41[,vars.G])
m41.my.se


#let's test the convex hull: from -0.32 to 0.29, for gov litig under 234 case
cfactmat41 <- matrix(data = NA, nrow = 101, ncol = 10)
k <- -0.5
for(i in 1:101){
	cfactmat41[i,] <- c(k,1,k,1,1,k,k,mean(data$percham),0,0)
	k <- k + 0.01
	}	
chullq <- whatif(data = m41.glm, cfact = cfactmat41)
chullq$in.hull

#let's test the convex hull: from -0.32 to 0.10
cfactmat41 <- matrix(data = NA, nrow = 101, ncol = 10)
k <- -0.5
for(i in 1:101){
	cfactmat41[i,] <- c(k,1,k,0,0,0,0,mean(data$percham),0,0)
	k <- k + 0.01
	}	
chullq <- whatif(data = m41.glm, cfact = cfactmat41)
chullq$in.hull


vars.m41.M <- c("normnetwobs")
## Relevant double interaction
vars.m41.D1 <- c("normwobsgov")
vars.m41.D2 <- c("normwobs177")

## Triple interaction
vars.m41.T <- c("normwobsgov177")

## Interaction vars
vars.m41.11 <- c(vars.m41.M, vars.m41.D1, vars.m41.D2, vars.m41.T)
vars.m41.10 <- c(vars.m41.M, vars.m41.D1)
vars.m41.01 <- c(vars.m41.M, vars.m41.D2)
vars.m41.00 <- c(vars.m41.M)

## Quantities of interest
zeta.m41.11 <- m41.glm$coefficients[vars.m41.M] + m41.glm$coefficients[vars.m41.D1] + m41.glm$coefficients[vars.m41.D2] + m41.glm$coefficients[vars.m41.T]
zeta.m41.10 <- m41.glm$coefficients[vars.m41.M] + m41.glm$coefficients[vars.m41.D1]
zeta.m41.01 <- m41.glm$coefficients[vars.m41.M] + m41.glm$coefficients[vars.m41.D2]
zeta.m41.00 <- m41.glm$coefficients[vars.m41.M]

#standard errors
vcov.m41 <- my.vcov(m41.glm, data.m41[,vars.G])



grad.m41.11.eta <- c(0, 1, 0, 1, 0, 0, 1, 1, 0, 0,0)
grad.m41.10.eta <- c(0, 1, 0, 1, 0, 0, 0, 0, 0, 0,0)
grad.m41.01.eta <- c(0, 1, 0, 0, 0, 0, 1, 0, 0, 0,0)
grad.m41.00.eta <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0,0)

se.m41.11 = sqrt(grad.m41.11.eta %*% vcov.m41 %*% grad.m41.11.eta)
se.m41.10 = sqrt(grad.m41.10.eta %*% vcov.m41 %*% grad.m41.10.eta)
se.m41.01 = sqrt(grad.m41.01.eta %*% vcov.m41 %*% grad.m41.01.eta)
se.m41.00 = sqrt(grad.m41.00.eta %*% vcov.m41 %*% grad.m41.00.eta)

## Put in a vector for convenience
t5m1 <- c(zeta.m41.11,se.m41.11,zeta.m41.10,se.m41.10,zeta.m41.01,se.m41.01,zeta.m41.00,se.m41.00)

